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C\ ■ Abstract 

5 : 

r- 1 . The QCD phase diagram may feature a critical end point at a 

Q_i| temperature T and baryon chemical potential /i which is accessible 

O 1 in heavy ion collisions. The universal long wavelength fluctuations 

which develop near this Ising critical point result in experimental 
signatures which can be used to find the critical point. The magnitude 
of the observed effects depends on how large the correlation length £ 
^ ■ becomes. Because the matter created in a heavy ion collision cools 

through the critical region of the phase diagram in a finite time, critical 
slowing down limits the growth of £, preventing it from staying in 
equilibrium. This is the fundamental nonequilibrium effect which must 
be calculated in order to make quantitative predictions for experiment. 
We use universal nonequilibrium dynamics and phenomenologically 
motivated values for the necessary nonuniversal quantities to estimate 
how much the growth of £ is slowed. 
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1 Introduction 



1.1 The Critical Point 

One goal of relativistic heavy ion collision experiments is to explore and map 
the QCD phase diagram as a function of temperature and baryon chemical 
potential. Recent theoretical developments suggest that a key qualitative 
feature, namely a critical point which in a sense defines the landscape to be 
mapped, may be within reach of discovery and analysis by the CERN SPS 
or by RHIC, if data is taken at several different energies |l|, [|. The discovery 
of the critical point would in a stroke transform the map of the QCD phase 
diagram from one based only on reasonable inference from universality, lattice 
gauge theory and models into one with a solid experimental basis [[|. 

In QCD with two massless quarks {m u ^ = 0; m s = oo) the phase 
transition at which chiral symmetry is restored is likely second order and 
belongs to the universality class of 0(4) spin models in three dimensions 
Below T c , chiral symmetry is broken and there are three massless pions. At 
T = T c , there are four massless degrees of freedom: the pions and the sigma. 
Above T = T c , the pion and sigma correlation lengths are degenerate and 
finite. 

In nature, the light quarks are not massless. Because of this explicit 
chiral symmetry breaking, the second order phase transition is replaced by 
an analytical crossover: physics changes dramatically but smoothly in the 
crossover region, and no correlation length diverges. This picture is consistent 
with present lattice simulations 0, which suggest T c ~ 140 — 170 MeV ||. 

Arguments based on a variety of models |7|, |S], |9], [TO], |TT], |TJ, [T3|, |TJ] indicate 



that the transition as a function of T is first order at large /i. This suggests 
that the phase diagram features a critical point E at which the line of first 
order phase transitions present for fi > fiE ends, as shown in Figure l.Q 
At he, the phase transition at T = T% is second order and is in the Ising 



universality class |Tj], [T2[ . Although the pions remain massive, the correlation 
length in the a channel diverges due to universal long wavelength fluctuations 
of the order parameter. This results in characteristic signatures, analogues 
of critical opalescence in the sense that they are unique to collisions which 
freeze out near the critical point, which can be used to discover E [0, §]. 
The position of the critical point is, of course, not universal. Furthermore, 



1 If the up and down quarks were massless, E would be a tricritical point, at which the 
first order transition becomes second order. 
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Figure 1: Sketch of the QCD phase diagram as a function of temperature T 
and baryon chemical potential //. Chiral symmetry is broken at low T and \x. 
As T is increased, chiral symmetry is approximately restored via a smooth 
crossover to the left of E or a first order phase transition to the right of E. 
The symmetry is only approximately restored because the light quarks are 
not massless. At the critical point E at which the line of first order phase 
transitions ends, the transition is second order and is in the Ising universality 
class. (At large fi and small T, there are color superconducting phases which 
we do not discuss in this paper.) The Ising model r-axis and /i-axis and the 
trajectories a, b and c will be discussed in Section 2. 



it is sensitive to the value of the strange quark mass. /j,e decreases as m s 
is decreased |]], and at some m c s , it reaches he = and the transition 
becomes entirely first order [15]. The value of m c s is an open question, but 
lattice simulations suggest that it is about half the physical strange quark 
mass JIB], |17], although these results are not yet conclusive JTj|. Of course, 
experimentalists cannot vary m s . They can, however, vary fi. The AGS, 
with beam energy 11 AGeV corresponding to yfs = 5 GeV, creates fireballs 
which freeze out near /i ~ 500 — 600 MeV E9|. When the SPS runs with 
^Js = 17 GeV (beam energy 158 AGeV), it creates fireballs which freeze out 
near /i ~ 200 MeV [O]. RHIC will make even smaller values of fi accessible. 



By dialing i/i and thus /i, experimenters can find the critical point E. 
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1.2 Detecting the Critical Point 

Predicting and thus suggesting the y/s to use to find E, is beyond 
the reach of present theoretical methods because /i£ is both nonuniversal 
and sensitively dependent on the mass of the strange quark. Crude models 
suggest that he could be ~ 600 — 800 MeV in the absence of the strange 
quark JCT , [12j] ; this in turn suggests that in nature /i e may have of order half 



this value, and may therefore be accessible at the SPS if the SPS runs with 
y/s < 17 GeV. However, at present theorists cannot predict the value of \xe 
even to within a factor of two. The SPS can search a significant fraction of 
the parameter space; if it does not find E, it will then be up to the RHIC 
experiments to map the < 200 MeV region. 

Locating E on the phase diagram can only be done convincingly by 
an experimental discovery. Theorists can, however, do reasonably well at 
describing the phenomena that occur near E, thus enabling experimenters 
to locate it. This is the goal of Ref. [0]. The signatures proposed there 
are based on the fact that E is a genuine thermodynamic singularity at 
which susceptibilities diverge and the order parameter fluctuates on long 
wavelengths. The resulting signatures are nonmonotonic as a function of y/s: 
as this control parameter is varied, we should see the signatures strengthen 
and then weaken again as the critical point is approached and then passed. 

The simplest observables to use are the event-by-event fluctuations of the 
mean transverse momentum of the charged particles in an event, px, and 
of the total charged multiplicity in an event, N. One analysis described 
in detail in Ref. || is based on the ratio of the width of the true event- 
by-event distribution of the mean p? to the width of the distribution in a 
sample of mixed events. This ratio was called y/F- NA49 has measured 
\/F = 1.002 ± 0.002 [^, H|, which is consistent with expectations for 
noncritical thermodynamic fluctuations.^ Critical fluctuations of the a field, 
i.e. the characteristic long wavelength fluctuations of the order parameter 
near E, influence pion momenta via the (large) amr coupling and increase 
y/~F p|. The effect is proportional to £f reezeout , where £f r eezeout is the a- field 

2 In an infinite system made of classical particles which is in thermal equilibrium, 
y/F = 1. Bose effects increase y/F by f — 2% |2f], ||; an anticorrelation introduced 
by energy conservation in a finite system — when one mode fluctuates up it is more likely 
for other modes to fluctuate down — decreases y/F by 1 — 2% Q; two-track resolution 
also decreases y/F by f — 2% [^0). The contributions due to correlations introduced by 
resonance decays and due to fluctuations in the flow velocity are each significantly smaller 
than 1% §. 
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correlation length of the long- wavelength fluctuations at freezeout 0. If 
£freezeout ~ 6 fm, the ratio ^/F increases by 10 — 20%, fifty times the statistical 
error in the present measurement [0. This observable is valuable because 
data on it has been analyzed and presented by NA49, and it can therefore be 
used to learn that Pb+Pb collisions at 158 AGeV do not freeze out near E. 

Once E is located, however, other observables which are more sensitive to 
critical effects will be more useful. For example, a i/^soft, defined using only 
the softest 10% of the pions in each event, will be much more sensitive to the 
critical long wavelength fluctuations. The higher pt pions are less affected 
by the a fluctuations 0, and these relatively unaffected pions dominate the 
mean of all the pions in the event. This is why the increase in \fF near 
the critical point will be much less than that of \/F so f t . 

The multiplicity of soft pions is an example of an observable which may 
be used to detect the critical fluctuations without an event-by-event analysis. 
The post-freezeout decay of sigmas, which are copious and light at freezeout 
near E and which decay subsequently when their mass increases above twice 
the pion mass, should result in a population of pions with px ~ which 
appears only for freezeout near the critical point 0. If £f ree zeout ^ l/ m T> this 
population of unusually low momentum pions will be comparable in number 
to that of the "direct" pions (i.e. those which were pions at freezeout) and 
will result in a large signature. 

The variety of observables which should all show nonmonotonic behavior 
near the critical point is sufficiently great that if it were to turn out that 
He < 200 MeV, making E inaccessible to the SPS, all four RHIC experiments 
could play a role in the study of the critical point. 

1.3 How Large Can £ Grow? 

Our purpose in this paper is to estimate how large £f ree zeout can become, thus 
making the predictions of Ref. [§] for the magnitude of various signatures 
more quantitative. In an ideal system of infinite size which was held at 
T = Te] fi = [j,e for an infinite time, the correlation length £ would be 
infinite. Ref. |2| estimated that finite size effects limit £ to be about 6 fm 
at most. We will argue in this paper that limitations imposed by the finite 
duration of a heavy ion collision are more severe, preventing £ from growing 
larger than about 2/Tg ~ 3 fm. 
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1-4 T E , T freezeout , and T 



We will do the calculation in the next section in the three-dimensional Ising 
model, as appropriate for describing the universal dynamics of the long 
wavelength fluctuations near the critical point. However, in order to relate a 
calculation in the Ising model to experiments which explore the QCD phase 
diagram, we will need numerical values for three temperature scales. Several 
other nonuniversal quantities will also enter our calculation; we will discuss 
them in the next section as they arise. We will see that in the end, only one 
combination of nonuniversal quantities plays a role in our estimates. 

We expect Te to be slightly less than the temperature range at which the 
crossover occurs at fi = 0. We therefore take T# = 140 MeV, at the low end 
of lattice estimates for the /i = crossover temperature. 

As we have discussed at length, we know very little about \ie- Fortunately, 
we will not need a numerical value for /i£ below. 

Pb+Pb collisions at 158 AGeV freeze out at about 120 MeV, and NA49 
data [2(| demonstrate clearly that they do not freeze out near E ||. We also 
know [I| that if the matter produced in a heavy ion collision comes near E, 
the large specific heat characteristic of E will cause the system to "linger" - 
the expansion will cause the energy density to decrease as usual, but this will 
result in an unusually slow temperature decrease. The freezeout temperature 
is therefore expected to be unusually close to the critical temperature for 
collisions which have the appropriate // to pass near E. For concreteness, 
we will take Tf reezeout = 130 MeV.0 (If the freeze-out temperature in Pb+Pb 



collisions at 158 AGeV is closer to 100 MeV, as some authors estimate f2H 
then it may be better to estimate that collisions which pass near E freezeout 
at T freezeout = 120 MeV.) 

Finally, we need to estimate T , the temperature at which we can begin an 
Ising model treatment. The three dimensional Ising model is only valid close 
enough to E that the correlation length £ > 1/Te- In this critical region, 
the long wavelength fluctuations of the order parameter become effectively 
three dimensional. (We will find that £ is never 3> 1/Tg. This means that 
our estimates are not precise.) We need to know how far above Te the 



equilibrium correlation length is larger than 1/Tg. The model of Ref. flT 
suggests that £ eq > 1/T& for (T — Te)/Te < 0.2 — 0.4. This estimate is 
based on a mean field analysis of a toy model, and so should not be taken 
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Note that experimenters do have some control over Tfreezeout- Using smaller ions 



results in a fireball which freezes out earlier, at a larger Tfreezeout 22 
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too seriously. For concreteness we shall assume that £ eq = l/T E = £ at 
T = 180 MeV, 40 MeV above T E ~ 140. We will use £ = 1-4 fori to set the 
scale below, in the sense that we will estimate the factor by which £/£o grows 
as the system cools. £o = 1-4 fori is simply a definition; To, the temperature 
at which the equilibrium correlation length £ eq = £o> is a quantity which must 
be estimated and which will affect our results. 



2 Slowing Out of Equilibrium 

The nonequilibrium dynamics which we analyze in this paper is fundamental 
in the sense that it is guaranteed to occur in a heavy ion collision which passes 
near E, even if local thermal equilibrium is achieved at a higher temperature 
during the earlier evolution of the plasma created in the collision. We assume 
early thermal (although not necessarily chemical) equilibration, and ask how 
the system evolves out of equilibrium as it passes E. More precisely, we will 
assume that when the system has cooled to T = To = 180 MeV, it is in 
equilibrium, with £(T ) = £ eq (To) = £ - For the present, assume that the 
system cools through the critical point E, as sketched in trajectory (a) of 
Figure 1. If it were to cool arbitrarily slowly, £ = £ eq would be maintained 
at all temperatures, and £ would diverge at Tg. However, it would take 
an infinite time for £ to grow infinitely large. Indeed, near a critical point, 
the long correlation length results in long equilibration times, a phenomenon 
known as critical slowing down. This means that the correlation length 
cannot grow as fast as £ eq , and the system cannot stay in equilibrium. 

We describe the effects of critical slowing down on the time development 
of the correlation length using the following equation for m a (t) = 

s m ' w= - r(m ' (t)) ( m ' (t) "d(i)) ■ (2J) 

Here, F parametrizes the rate at which an out-of-equilibrium value of m a 
approaches its equilibrium value. If m a is close to its equilibrium value, the 



theory of dynamical critical phenomena [25| tells us that 



r(m CT ) = ^ (m^ ) z (2.2) 

where z is a universal exponent and we have used £o to set the scale, making 
A a dimensionless constant. Knowing that we are interested in a system 
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which is in the same static universality class as the 3-dimensional Ising 
model is not enough to tell us z. There are in general several different 
dynamical universality classes corresponding to a given static universality 
class. However, knowing in addition that: (i) the chiral order parameter 
is not a conserved quantity; (ii) there are other conserved quantities in the 
system, like the baryon number density; and (iii) there are no Poisson bracket 
relations between the order parameter and the conserved quantities, tells us 
that our system belongs in the dynamical universality class named Model C in 
Halperin and Hohenberg's classification |2f| of dynamical critical phenomena, 
and has 

z = 2 + a/i/^2.17, (2.3) 



where we have taken a = 0.11 and v = 0.630 from Ref. p6[ . The 
dimensionless constant A is nonuniversal. We have no way to estimate it 
other than to guess that it is of order 1. We will explore the sensitivity of 
our results to different choices of A below. 

We will use the differential equation ( |2.1| ) to analyze how critical slowing 
down prevents the correlation length £ from "tracking" £ eq (T(t)). Critical 
slowing down guarantees that the system falls out of equilibrium. Note that 
the differential equation has only been derived for small departures from 
equilibrium; once rn a —^~,} is not small, its use is not quantitatively justified.^ 

We have initial conditions for the differential equation (|2 . 1|) , namely 
m a (0) = l/£o- Therefore, all we need in order to solve it is a description 
of £ eq (£)- This requires £ eq (T), which we discuss below, and also requires a 
description of the cooling T(t). This can be estimated using hydrodynamic 
and cascade model calculations, although these describe T(t) assuming the 
plasma is not cooling near the critical point E. Hydrodynamic models (see, 
e.g., Refs. p|, pTfl ) describe T{t) at central rapidity in the center of mass 



frame via 

£ = --■ (2-4) 
dt Kt v 1 

Since we are only interested in a relatively small range of temperatures around 
T E , it will suffice for us to treat dT/dt as constant in time. We discuss the 

4 For example, one might try the equation d^/dt = — r(£ — £ q) ; instead of the equation 



(2.1) for dm a /dt. These two equations give the same results for small departures from 
equilibrium, but they do not agree in all circumstances. For example, in a system which is 
not cooling and which has T — Te and £ cq = oo for all time, only (2.1) yields the correct 
result, namely 77i CT (<) ~ t~ x l z at late time. 
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effects of the time dependence of dT / dt below. The expression ( |2.4j ) assumes 
that the T-dependence of the energy density is e ~ T K as in a resonance gas, 
for which k ~ 6 ||28|| . The timescale to is not constant over the whole history of 
the collision. A simplified estimate (made by equating to with the scattering 
time) suggests that in Pb+Pb collisions at 158 AGeV, to is between 4 and 10 
fm at times of interest to us |^7||. This suggests —dT/dt = (2 — 6) MeV/fm 
at T = 140 MeV. Careful analysis favors to closer to 4 fm |[24!| . This agrees 
with a recent analysis of these collisions using the URQMD cascade model, 
which suggests -dT/dt « 5 MeV/fm at T = 140 MeV p|. These estimates 
are all for cooling through T = 140 MeV at a \i such that one is not near the 
critical point. As we discussed above, the cooling rate is likely to be unusually 
low near E because of the large specific heat there; we will therefore take 
—dT/dt ~ 4 MeV/fm as our estimate, noting also that the cooling rate at 
RHIC will be slower still. 

We wish to use the three-dimensional Ising model to describe £ e q(^, /x) 
near E. In the Ising model, the order parameter M (the magnetization) 
and the correlation length £ are functions of the reduced temperature r and 
the magnetic field h. (In the Ising model, r is defined as (T — T c )/T c and 
is usually called t; we reserve the symbol t for time, however.) The critical 
point is at r = h = 0; at this point, M = and £ eq = oo. For r < 0, there is 
a first order phase transition as a function of h at h = between M = \r\@ 
for h = 0+ and M = —\r\@ for h = 0—. The exponent is (3 = 0.326 for the 
3-dimensional Ising model J2EJ. For r > 0, M increases smoothly through 
zero as h goes from negative to positive. For r = 0, the order parameter is 
M = sign(/i)|/i| 1 /' 5 , with 5 = 4.80 |6 



We can now discuss how the Ising model r- and h-ax.es are mapped onto 
the (T, jj) plane. The r-axis is the direction tangential to the line of first 
order phase transitions ending at E. This is shown in Figure 1. There is 
no guarantee that the Tt-axis is perpendicular to the r-axis when both are 
mapped onto the (T, //) plane. This mapping will in general deform the Ising 
axes, but we have no way of estimating this deformation.^ For simplicity, we 

J In the clectroweak phase diagram as a function of T and Higgs mass to#, there is 
also a line of first order phase transitions ending at an Ising critical point. Here, the 
explicit mapping between Ising axes and the (T, m#) plane has been constructed 
This is possible only because there are reliable numerical methods for analyzing the full, 
nonuniversal theory in the (T, mii) plane. Universality arguments alone, which is all that 
we have at our disposal in the absence of lattice simulations at nonzero do not tell us 
how the Ising axes should be deformed in the (T, fi) plane of Figure 1. 
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draw the /i-axis perpendicular to the r-axis in Figure 1. In thinking through 
the mapping between QCD and the Ising model even qualitatively, it is 
important to note that the QCD order parameter ( the chiral condensate (qq)) 
is offset with respect to the Ising model order parameter (the magnetization 
M). In the Ising model, M = at the critical point and along the first 
order line one has phase coexistence between phases whose M's are equal in 
magnitude and opposite in sign. In QCD, (qq) ^ at the critical point E, 
because of the explicit breaking of the 0(4) symmetry by quark mass terms. 
Near E, (qq) corresponds to M plus an offset, and the phase coexistence is 
between phases with differing values of (qq) which both have the same sign. 
In Figure 1, we take the —h side of the Ising coexistence line to correspond to 
the higher temperature side of the QCD coexistence line, so that increasing 
M corresponds to increasing the magnitude of (qq). 

The matter created in heavy ion collisions at SPS energies will follow a 
trajectory in the (T, /i) plane which is approximately vertical as it cools. (See, 
for example, Ref. ||29|| .) We therefore begin by considering trajectory (a) of 
Figure 1, which follows the h this is likely not a bad approximation 

to cooling at almost constant We have analyzed trajectories which pass 
through E at a variety of angles, for example like trajectory (b) in Figure 1. 
The results do not differ qualitatively from those we present in detail for a 
trajectory along the /i-axis, unless the trajectory passes through E almost 
parallel to the r-axis. At the end of this section, we will present results for 
trajectories like (c) in Figure 1, which miss E but come close to it. 

Let us take the initial temperature in our calculation, T = To = 180 MeV, 
to correspond to h = ho = —0.2. Along the /i-axis,[| as in trajectory (a), the 
equilibrium correlation length is a power law in h: 

-u/138 

(2.5) 

where we have normalized £ by setting £ e q(^o) = 1- That is, we measure £ 
in units of £o = 1-4 fm. With units chosen, we can now rewrite the equation 

6 Note that the r-direction, corresponding to the reduced temperature direction in the 
Ising model, is almost perpendicular to the T direction in QCD. This is another reason 
why we have labeled it by a letter other than t. 

'Along the r-axis, £ C q ~ r~ v \ along the /i-axis, £ oq ~ h~ u ^ s ; for trajectories like (b) 
of Figure 1 which pass through E at generic angles, the larger exponent (v/(35) is the 
relevant one. Our /i-axis analysis is therefore a good guide to the generic case. 



h 
h 



9 



(2~T) which describes the dynamics of the growth of the correlation length 



£ = l/m a in terms of Ising model variables as 
d 



dh mAh) = -a^(h)) (^)-^)J (2-6) 

where the nonuniversal constant a is related to the other nonuniversal 
parameters we have discussed by: 

a = A (fY 1 = A(-!a-^)' 1 . (2.7) 

\dtl \T -T E dt y ' 



Nonuniversal parameters appear in equation (|2.6f) only in the single 
combination a. Taking the nonuniversal constant from ( |2.2j ) to be A ~ 1, 
using (T - T E ) = 40 MeV, dT/dt = -4 MeV/fm and h = -0.2 yields the 
estimate 

a ~ 50 . (2.8) 

In fact, because £ e q(^) is a power law in h, if one changes h$ and then redefines 
the units of £ so that £ eq (^o) is again set to one, equation (|2.6| ) is unaffected. 
Our results are therefore determined solely by A, (To — Te) and dT/dt in the 
single combination a, together with the assumption that the system begins 
in equilibrium at T = T . 

We can now use ( |2.6|) to learn how much £ grows relative to £o — 1-4 fm. 
Given the uncertainties in the determination of a, in Figure 2 we show £(/i) 



obtained by solving ( |2.6| ) for a = 25, 50, 100. Four lessons are apparent: 

First, critical slowing down has a large effect. Although by assumption 
we begin in thermal equilibrium with £ = £ eq at T = To, the fact that the 
dynamics slows down in the vicinity of E prevents £ from tracking £ eq and 
growing very large. 

Second, our results do not depend sensitively on the parameter a. This is 
fortunate, since there are so many uncertainties involved in estimating a. For 
a = 25,50, 100, the maximum correlation length which is achieved is 1.8£o, 
2.1£ , 2.5£ , corresponding to 2.6, 3.0, 3.4 fm. This means that although our 
estimate is only qualitative, it is clear that £ cannot grow as large as 6 fm. 
(To obtain a maximum value of £ = 4£ would require a = 1000. Although a 
is uncertain, this large a value seems out of the question.) We estimate that 
£ grows to about twice £o ; corresponding to approximately 3 fm. 

Third, since previous work || suggests that finite size effects limit £ to 
£ < 6 fm, we conclude that slowing out of equilibrium (i.e. the combination of 
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-0.2 -0.15 -0.1 -0.05 0.05 0.1 

Figure 2: Behavior of the correlation length for cooling through the critical 
point along the fa-axis of Figure 1. The equilibrium correlation length is 
shown as a dashed line. The true correlation length is shown for (bottom 
to top) a = 25, 50, 100. Our units, described in the text, are such that 
h = -0.2,-0.1,0,0.1 corresponds to T = 180,160,140,120 MeV, and f is 
measured in units of £ — 1-4 fin. 

finite time and critical slowing down) imposes the more stringent constraint 
on £. We have analyzed this nonequilibrium effect as if the system were 
spatially homogeneous. Had we found correlation lengths growing beyond 
6 fm, we would have to do a much more complicated analysis, taking both 
the finite time and the inhomogeneous spatial dynamics into account. Since 
we find that £ only grows to about 3 fm, this is not necessary. 

Fourth, just as critical slowing down prevents £ from growing as fast as 
£ cq does, it also prevents £ from shrinking as fast as £ cq does after E has 
been passed. Whether we estimate Tf reez eout ~ 130 MeV (corresponding to 
h = .05) or ~ 120 MeV (h = 0.1) for trajectories passing near E, one finds 
tfreezeout ~ 2£o- We can also argue that even if an a as large as 1000 were 
possible, our conclusions would be little affected: If a = 1000, and £ follows 
£ cq closely enough that it increases to 4£ , £ also tracks £ eq more closely as 
it decreases below T E . If a = 1000, it turns out that £ is quite close to £ eq 
by the time h = .05. Thus, although increasing a to a ridiculous extent does 
increase the maximum value of £, it has little effect on £f ree zeout- This is further 
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evidence that although our estimate £f re ezeout ~ 2£ ~ 3 fm is qualitative, it 
is robust. 

The a-dependence of the maximum value attained by £ can be understood 
analytically at large a. For large a, m a tracks its equilibrium value m£. q = 
\h/ho\ v ^ s well until \h/ho\ is quite small. If we define e = m a — m^. q , we can 
use ( j2.6| ) to show that e < m£ q as long as 

h ( v 1 \ i+W^ 

> -sir • ( 2 - 9 ) 



/io \(38 a 

If we assume that once £ begins to drop out of equilibrium (i.e. once e begins 
to grow comparable to m^. q ) , little further growth of £ occurs before £ reaches 
its maximum, we predict that £ will peak at 

£ max = ( C a)^W = ( ca ) - 215 , (2.10) 



for some constant c. The maxima of solutions to ( j2.6| ) obtained numerically 
follow this scaling relation (with c = 0.65) quite accurately once a > 1000 or 
so. Even at much smaller a, as in Figure 2, £ max is within a few percent 
of that in ( [2 .ICQ . This scaling relation explains why our results are so 
weakly dependent on a. Note that even with the scaling relation in hand, 
full solutions as in Figure 2 are of value because they allow us to estimate 

ffreezeout and not JUSt £ 

max ■ 

We have to this point assumed that dT/dt is approximately constant 
as the system cools through Te- This is an oversimplification. It is more 
reasonable to assume that ds/dt is approximately constant, where s is the 
entropy density. Since 

dT _ 1 ds 

~dJ~c^dt ^" ll > 

and the specific heat Cy is peaked at Te, we expect that dT/dt is unusually 
small near Te- As we discussed above, this "lingering" results in a Tf reez eout 
which is unusually close to T E [0. Here, we estimate the effect of lingering 
near E on the growth of £. Along the h-axis, the specific heat due to the long 
wavelength sigma fluctuations diverges like Cy ~ h" 1 ^ 5 ~ £2q in thermal 
equilibrium Q. The exponent 7 = 1.240 |2f|. We take Cy — c\ + 02^^, 



where c± is the specific heat due to all the degrees of freedom other than the 
sigma and is smooth near Te- Note that Cy depends on the actual correlation 
length £, and not on £ eq . In our dynamical nonequilibrium setting, therefore, 
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Cy peaks but does not diverge. We can implement lingering in our calculation 
by replacing the constant a in (|2.6| ) by 



a(h) =a (1-b) +b[m (7 (h)}-~ t/ " . (2.12) 

Here, a is the same constant as before and the constant b = €2/(01 + c^) is 
the fraction of the specific heat at T = T which is due to sigma fluctuations. 
This fraction is perhaps about 0.1 and is surely less than 0.25. As the system 
cools from T to Te, the sigma contribution to Cy grows and peaks. We find 
that changing constant a to a(h) as in ( |2.12| ) with b = 0.25 increases £ n 



by about 10% beyond that shown in Figure 2, and increases £f ree zeout by 
somewhat less. For b = 0.1, the increase in £ max is about 5%. We conclude 
that because Cy receives contributions from all degrees of freedom and not 
just from the sigma fluctuations, and because Cy, like £, peaks but does not 
diverge, the reduction in dT/dt near E is not large enough to significantly 
increase £ beyond our previous estimates. 

We now ask how much our results change if we consider trajectories like 
(c) in Figure 1 which come close to, but miss, E. Our analysis can easily 
be extended to cover those trajectories which pass E on the crossover side 
(r > 0; T < Te). In an appendix, we present the Ising model expression for 
Ceq( r ,h) near the critical point. We use this expression to evaluate for 
trajectories parallel to the /i-axis with r = 0.12, r = 0.19 and r = 0.33, for 
which £ eq peaks at 4£ , 3£ and 2£ . The results are shown in Figure 3 in 
which we have taken a = 50. Note that in plotting Figure 3 we have defined 
£0 = ho) — 1 anew for each r. We see that as long as £ eq peaks at 3£o or 
higher, the dynamics of £ is almost the same as for the trajectory of Figure 
2 which goes precisely through E. Even for a trajectory which misses E by 
enough that £ eq peaks at only 2£ , the actual correlation length £ grows by 
a factor which is within 20% of that for trajectories which pass arbitrarily 
close to E. Just as the growth of £ is robust with respect to changes in a, 
it is robust with respect to how close the trajectory comes to E, for those 
trajectories which come close enough. 

Our analysis is not sufficient to describe the dynamics for those 
trajectories which pass to the first order side of E, because we do not treat 
the dynamics of bubble nucleation and phase coexistence. Near enough to E, 
though, the first order transition is so weak that it will not have detectable 
effects given the finite length and time scales in a heavy ion collision, and 
the physics is likely qualitatively similar to that we have analyzed on the 
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Figure 3: The dashed curves show £ eq (/i) for trajectories with (top to bottom) 
r = 0.12, 0.19, 0.33. The solid curves show the corresponding nonequilibrium 
correlation lengths £, assuming a = 50. Curves with different r have each 
been normalized to begin at £,(h ) = 1. 

crossover side of E. Farther from E on the first order side, this is not the 
case. Farther from E, though, the correlation length is never large. 

3 Consequences 

Our results have a number of consequences which should be taken into 
account both in planning experimental searches for the QCD critical point, 
and in planning future theoretical work. 

Because of critical slowing down, the correlation length in a heavy ion 
collision cannot grow as fast as it would in equilibrium; this means that 
Cfreezeout is likely about 3 fm for trajectories passing near E. Although 
finite size effects alone would allow a correlation length as large as 6 fm, 
this is unrealistic to expect in a heavy ion collision. This effect arises due 
to guaranteed nonequilibrium physics: even if heavy ion collisions achieve 
local thermal equilibrium above the transition, as we have assumed, if they 
cool through the transition near the critical point they must "slow out of 
equilibrium." By this we mean that the correlation length cannot grow as it 
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would in equilibrium, because the long wavelength dynamics are slow near E. 

Critical slowing down also prevents the correlation length from decreasing 
quickly after passing the critical point. One therefore need not worry about 
£ decreasing significantly between the phase transition and freezeout. 

One need not hit E precisely in order to find it. The results shown in 
Figure 3 demonstrate that if one were to do a scan with collisions at many 
finely spaced values of the energy and thus /x, one would see signatures of 
E with approximately the same magnitude over a broad range of /x. The 
magnitude of the signatures will not be narrowly peaked as /x is varied. As 
long as one gets close enough to E that the equilibrium correlation length 
is (2 — 3)£o5 the actual correlation length £ will grow to ~ 2£ - There is no 
advantage to getting closer to E, because critical slowing down prevents £ 
from getting much larger even if £ eq does. Data at many finely spaced values 
of /x is not called for.f] 

Only one combination of the nonuniversal quantities (called a above) 
plays an important role in estimating the dynamics of £. The uncertainty 
in a is the sum of that in its three factors: A (the nonuniversal constant in 
the dynamical scaling law fl2.2|) ), dT/dt and T — T#. It is already fortunate 
that only one combination a matters; it is even more fortunate that our 
results are not very sensitive to the value of a. This means that although our 
results are not completely quantitative, they are robust. In addition to the 
uncertainty in a, however, our results cannot be treated as precise because 
the QCD dynamics are precisely described by the three-dimensional Ising 
model dynamics only if £ ^> 1/Tg, and we have found that £ does not grow 
beyond ~ 2/Tg. 

There are a number of steps that could be taken in future work to refine 
our estimate. One could do a more complete job of analyzing the universal 
dynamics of a system which passes near an Ising critical point. For example, 
instead of simply writing a differential equation for £, one could follow the 
full 3+1-dimensional dynamics in a Langevin simulation, from which one 
would measure £. Doing this, however, would still leave one facing the same 
nonuniversal uncertainties which we face in our treatment. If we simply 

8 Analysis within the toy model of Ref. |Q suggests that in the absence of the strange 
quark, the range of fi over which £ oq > 2 fm is about A/x ~ 120 MeV for he ~ 800 MeV. 
Similar results can be obtained pl| within a random matrix model |l2] | . It is likely over- 
optimistic to estimate A/x ~ 120 MeV when the effects of the strange quark are included 
and he itself is reduced. A conservative estimate would be to use the models to estimate 
that Ah/he ~ 15%. 
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ask how to reduce the uncertainty in a, perhaps the hardest part of this 
task would be a reliable calculation of A, as that would require a reliable 
calculational method for QCD dynamics at nonzero T and /i. The other two 
ingredients in a are likely to become better known as the modeling of heavy 
ion collisions and the analysis of data from these collisions proceed. It seems, 
though, that the uncertainty in A will prevent a fully quantitative calculation 
of a for the foreseeable future. Our results are sufficiently insensitive to a that 
they suffice to estimate the magnitude of signatures; when these signatures 
are found, perhaps they will give us more quantitative information about the 
nonuniversal quantities which go into a. 

Knowing that we are looking for £f rce zcout ~ 3 fm is very helpful in 
suggesting how to employ the signatures described in detail in Ref. 0. The 
excess of pions with ~ arising from post-freezeout decay of sigmas 

is large as long as £f re ezeout ~ l/m-m and does not increase much further 
if £freezeout is longer. This makes it an ideal signature. The increase in the 
event-by-event fluctuations in the mean transverse momentum of the charged 
pions in an event (described by the ratio y/~F of Ref. [0]) is proportional to 
£freezeout- The results of Ref. @ suggest that for £f ree zeout ~ 3 fm, this will be 
a 3 — 5% effect. This is ten to twenty times larger than the statistical error 
in the present NA49 data, but not so large as to make one confident of using 
this alone as a signature for E. The solution is to use signatures which focus 
on the event-by-event fluctuations of only the low momentum pions. Unusual 
event-by-event fluctuations in the pion momenta arise via the coupling 
between the pions and the sigma order parameter which, at freezeout, is 
fluctuating with correlation length £f ree zeout- This interaction has the largest 
effect on the softest pions 0. v^soft, described in the introduction, is a 
good example of an observable which takes advantage of this. Depending 
on the details of the cuts used to define it, it should be enhanced by many 
tens of percent in collisions passing near E. Ref. @ suggests other such 
observables, and more can surely be found. Together, the excess multiplicity 
at low momentum (due to post-freezeout sigma decays) and the excess event- 
by-event fluctuation of the momenta of the low momentum pions (due to their 
coupling to the order parameter which is fluctuating with correlation length 
^freezeout) should allow a convincing detection of the critical point E. Both 
should behave nonmonotonically as the collision energy, and hence /i, are 
varied. Both should peak for those heavy ion collisions which freeze out near 

E, With £ fre ezcout ~ 3 fm. 
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Appendix A Equilibrium Correlation Length 



In this appendix we present the equations that we used to calculate the 
equilibrium correlation length £ eq in the critical region as a function of 
reduced temperature r and external magnetic field h. We use Widom's 
scaling form |B2| 



e eq (r,M)=fM- 2 ^^ 



\M\ 



VP 



(A.l) 



in which M is the magnetization and v = 0.630 and (3 = 0.326 are the three- 
dimensional Ising model critical exponents [2B|. The e-expansion of g(x) is 



given in |32 to order e . 



g(x) = g £ {x) 



where z 



Q- 2v z{ 1 



36 



[(5 + 6 In 3)s -6(1 +2) In z) 



l + 2z 2 
72 



ln 2 z+(^-i)(l-ln3)- 



47 



56, 



(16z z z ) ln^ + 

216 V 3 3 7 

— f — + ^ 61 2 3 + 41 3 _ U) } 2 

216V 6 3 

-Lf61n"3 + ^ln3 + H + i /): 
216 V 3 9 3 



1 + 2 

^3 



and I 



1 In [g(l - x)} 
1 — x(l — x) 



(A.2) 



dx 



-2.344. 



/ in ( |A.1|) is a non-universal normalization constant, often set to one. Our 
choice of / and thus of units for £ is described in Section 2. The expression 
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Figure 4: Order parameter (magnetization) in the 3D Ising model as a 
function of reduced temperature r and applied field h. 



(A. 2) is valid everywhere on the (r, M) plane except the region of large x (or, 
equivalently, r ^> IM) 1 /^). The correct result at large x is 

^iargc(x) = > ( A - 3 ) 

and we therefore construct a function g(x) which smoothly interpolates 
between gwge^) at large x and g £ (x) at smaller x. The only remaining 
difficulty is at r > 0, M = 0. Although the scaling form (|A. 1|) with (|A.3|) is 
well-behaved in the M — >■ limit, and yields 



&q(r>0,M-0) = /|r| 



-V 

1 



at M = the scaling form is indeterminate and one must impose the 
condition £ eq (r > 0, M = 0) = ^ eq (r > 0, M -> 0). 

We want £, eq (r,h). With ^ eq (r, M) in hand, we must now obtain the 
magnetization M(r,h). The most convenient form for our purposes is the 
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Figure 5: Equilibrium correlation length as a function of reduced temperature 
r and applied field h. 



parametric equation of state (see [^6|, |32|, |33|]): 

M = M RP6 

h = h R^ 5 h(6) = h R ps (6 - 0.76201 6 3 + 0.00804 6 5 ) (A.4) 
r = R(l-6 2 ) 

Here M, r and h are parametrized in terms of the "radius" R > and "polar 
angle" 6.6 = corresponds to r > 0, h = 0; = ±1 corresponds to r = 
with positive and negative h respectively. h(6) is zero at 6 = ±1.154, which 
corresponds to r < 0, h = ±0. The function h(6) in ( |A.4| ) is from Ref. p6 



where it is constructed to be consistent with all that is known from the e- 
expansion, from perturbation theory, and from resummations thereof. We 
choose the normalization constants Mq and ho so that Mir = —1, h = 0) = 1 
and Mir = 0, h = 1) = 1. This guarantees that along the negative r axis 
M(r,h = 0) = \r\~P and along the h axis M{r = 0,h) = sign(h)\h\~ s . 
Numerically solving the last two equations of (|A.4j ) for R and 6 in terms of 
r and h, we use the first one to compute M(r, h), shown in Figure 4. This 
allows us to obtain £ eq (r, h), shown in Figure 5. 
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